package scu.maqiang.bes;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.DirectDSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class PossionCylinder {

	public static void main(String[] args) {
		Mesh2T3 mesh = new Mesh2T3();
		mesh.ring(0.5, 1, 256, 32);
		//mesh.displayBoundarys();
		//mesh.displayNodes();
		//mesh.displayBoundarys();
		BES2DT3L20 bs = new BES2DT3L20(mesh);
		//bs.displayMiddlePoint();
		bs.setDirichletBoundaryValue(5.0, 1);
		bs.setDirichletBoundaryValue(10.0, 2);
		double[][] H = new double[bs.GetNdof()][bs.GetNdof()];
		double[][] G = new double[bs.GetNdof()][bs.GetNdof()];
		bs.assembleLaplaceMatrix(new double[] {1.0}, H, G);
		double[] knowValue = bs.getKnownValue();
		double[] RHS = MVO.matmul(G, knowValue);
		//System.out.println(MVO.toString(knowValue));
		DirectDSolver solver = new DirectDSolver(H);
		double[] x = new double[bs.GetNdof()];
		solver.PLU();
		solver.PLUSolve(RHS, x);
//		System.out.println(Arrays.toString(x));
		double[][] result = bs.getSolution(x);
		//System.out.println(MVO.toString(result));
		//double[][] result2 = bs.computeDomainValues(result[0], result[1]);
		System.out.println("Result from boundary element method: ");
		System.out.println(bs.computeValue(result[0], result[1], null, new double[] {0.6, 0.0}));
		System.out.println(bs.computeValue(result[0], result[1], null, new double[] {0.7, 0.0}));
		System.out.println(bs.computeValue(result[0], result[1], null, new double[] {0.8, 0.0}));
		System.out.println(bs.computeValue(result[0], result[1], null, new double[] {0.9, 0.0}));
		double[][] domainResult = bs.computeDomainValues(result[0], result[1], null);
		mesh.toTecplot("BEMPossionResult.dat", domainResult);
		
		FES2T31 fs = new FES2T31(mesh);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		RHS = new double[fs.GetNdof()];
		x = new double[fs.GetNdof()];
		fs.assembleStiff(new double[] {1.0}, BVPType.COMMON, A);
		fs.applyBC_MBN(A, Direct.All, 1);
		fs.applyBC_MBN(A, Direct.All, 2);
		fs.applyBC_MBN(RHS, Direct.All, 5.0, 1);
		fs.applyBC_MBN(RHS, Direct.All, 10.0, 2);
		NewIterSSolver ssolver = new NewIterSSolver(A);
		ssolver.PCGSSOR(RHS, x, 1.5, 1);
		System.out.println("Result from finite element method: ");
		System.out.println(fs.getValue(x, new double[] {0.6, 0.0}));
		System.out.println(fs.getValue(x, new double[] {0.7, 0.0}));
		System.out.println(fs.getValue(x, new double[] {0.8, 0.0}));
		System.out.println(fs.getValue(x, new double[] {0.9, 0.0}));
	}

}
